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A time-domain numerical code based on the constitutive relations of nonlinear acoustics for simulating ultra¬ 
sound propagation is presented. To model frequency power law attenuation, such as observed in biological 
media, multiple relaxation processes are included and relaxation parameters are fitted to both exact fre¬ 
quency power law attenuation and empirically measured attenuation of a variety of tissues that does not 
fit an exact power law. A computational technique based on artificial relaxation is included to correct the 
non-negligible numerical dispersion of the numerical method and to improve stability when shock waves are 
present. This technique avoids the use of high order finite difference schemes, leading to fast calculations. 
The numerical code is especially suitable to study high intensity and focused axisymmetric acoustic beams 
in tissue-like medium, as it is based on the full constitutive relations that overcomes the limitations of the 
parabolic approximations, while some specific effects not contemplated by the Westervelt equation can be 
also studied. The accuracy of the method is discussed by comparing the proposed simulation solutions to 
one-dimensional analytical ones, to fc-space numerical solutions and also to experimental data from a focused 
beam propagating in a frequency power law attenuation media. 

PACS numbers: 43.58.Ta, 43.80.Sh, 43.35.Wa, 43.35.Fj, 43.25.Ts 


I. INTRODUCTION 

The significant development of ultrasound technology 
in the medical field in recent years has led to the need 
for simulation tools increasingly realistic. Effects like 
absorption in biological media, nonlinear propagation, 
heterogeneities, strong focusing, streaming, resonances, 
multiple scattering or the presence of discontinuities due 
to tissue layers or rigid boundaries have to be taken into 
consideration. The most general approach for ultrasound 
simulation is to directly solve the constitutive relations of 
the nonlinear acoustics. It also allows the explicit calcu¬ 
lation of the particle velocity, what can be used to com¬ 
pute important magnitudes as the vector components of 
the nonlinear acoustic intensity or the acoustic radiation 
force. 

The numerical resolution of the nonlinear constitutive 
equations in tissue-like medium supposes a difficult prob¬ 
lem due to the large size of the region of interest in 
relation to the size of the acoustic wavelength and the 
complexity of the model. Simplifying assumptions have 
been needed in the past for modeling beam patterns from 
ultrasound transducers, as one-way parabolic approx¬ 
imations, most based on the Khokhlov-Zabolotskaya- 
Kuznetsov equation (KZK) 1-8 . To overcome the valid¬ 
ity limitations of the parabolic approximations, i.e. for 
large aperture focused sound sources or modeling sound 
field near the acoustic source, many one-way numerical 
methods have been proposed, including phenomenologi¬ 
cal approaches 9 ’ 10 with tissue attenuation 11 , or based on 
the one-way formulation of the Westervelt equation 12,13 . 


Tissue inhomogeneity can be modeled in these one¬ 
way models 6 , like transmission though tissue layers with 
refraction, but they do not take into account backscat- 
tering and multiple reflections. More realistic models, 
e.g. those accounting for scattering from internal tissue 
structures, are based on the Westervelt-type full-wave 
equations 14-20 . This full-wave equation has been vali¬ 
dated for strongly focused sources 21 . However, due to 
the assumptions taken in the derivation of the Westervelt 
equation, the accuracy of this model is limited in prac¬ 
tical situations as (*) the modeling of rigid boundaries 
where the thermo-viscous boundary layer effects are not- 
negligible, i.e. in general case where the particle velocity 
field becomes rotational, (ii) situations where the second 
order Lagrangian density of acoustical energy not vanish, 
i.e. near the source or in situations where plane progres¬ 
sive waves does not exist and the acoustic held becomes 
complex due to multiple scattering, reverberation or res¬ 
onances, (in) situations where the equilibrium-state par¬ 
ticle velocity is not null, including the self generation of 
acoustic streaming. See Ref. Hamilton and Blackstock 22 
Chap. 3 for further discussion. 

The recent development of computational capacity 
has made possible to consider the full constitutive re¬ 
lations (i.e. without the assumptions discussed above). 
Thus, for small-amplitude acoustic waves, the linearized 
pressure-velocity formulation of constitutive relations in 
inhomogeneous media was solved by means of Finite- 
Differences in Time-Domain (FDTD) methods with fre¬ 
quency independent losses by Manry et al. 23 , or using 
two-step MacCormack finite-differences scheme by Mast 
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et al. 24 . Also, relaxation processes can be included in 
finite difference methods in an efficient way in order to 
model tissue attenuation and dispersion 25 . Furthermore, 
fc-space numerical methods have also been applied to 
solve the linearized first order equations in lossless in¬ 
homogeneous media 26 . In order to account for soft tissue 
losses, the computational solution of the fractional Lapla- 
cian by fc-space spectral methods have demonstrated to 
be extremely efficient due to the spatial frequency do¬ 
main representation of the acoustic field 2 '. 

In the case of nonlinear constitutive relations mod¬ 
els, the evolution of the acoustic magnitudes have been 
simulated in time-domain by means of finite differences 
schemes such as Dispersion Relation Preserving method 
(DRP) in ideal fluids and axisynnnetric domains 28 . 
Thermo-viscous losses in finite-differences methods have 
been widely used, see Sparrow et al. 29 . In order to in¬ 
troduce tissue attenuation in the governing equations 
time-dependent fractional derivatives can be included by 
convolutional operators. Thus, in Ref. Liebler et al. 30 
an efficient method has been presented, but although 
the memory requirements can be strongly reduced com¬ 
pared to direct convolutions the algorithm employs up 
to ten auxiliary fields and a memory buffer of three time 
steps. Furthermore, construction of specific causal mem¬ 
ory functions that models soft tissue attenuation and dis¬ 
persion in Navier-Stokes equations is also possible 31 , but 
certain time history must be stored in memory and in 
this case the computational domain was restricted to one 
dimensional propagation. 

In order to overcome those numerical limitations, 
recently fc-space and pseudo-spectral numerical meth¬ 
ods have been applied to constitutive relations in non¬ 
linear regime to solve fractional Laplacian operators 
efficiently 32 . Furthermore, in the case of domains of hun¬ 
dreds of wavelengths, when the cumulative phase error 
due to numerical dispersion of standard finite-difference 
schemes can not be neglected, those spectral numeri¬ 
cal methods have reported an improvement in accuracy 
of the numerical solution. This two factors, i.e. the 
negligible numerical dispersion and the efficient resolu¬ 
tion of fractional Laplacian operators, have led the spec¬ 
tral methods to be widely used in practical applications. 
However, their main limitation is that the implementa¬ 
tion of natural space discontinuities due to tissue layers 
or rigid boundary conditions leads to errors in the re¬ 
construction of the spectral information due to the poor 
convergence of Fourier series at jumps, i.e. the well- 
known Gibbs oscillations. Preventing this kind of errors 
is typically achieved by filtering the spatial spectrum 20 , 
so the theoretical spatial minimum sampling of two point 
per wavelength becomes larger. In addition, these errors 
propagate globally and affect to the accuracy all over 
the domain, in contrast with locally propagating errors 
in finite differences methods. On the other hand, taking 
into account the spatial discontinuity due to symmetry 
boundary condition, axisymmetric domains becomes not 
feasible by standard fc-space methods, and full 3D do¬ 


mains must be employed even for axisymmetric configu¬ 
rations. Those errors can be prevented by means of the 
recently developed Fourier Continuation (FC) method 33 . 
However, the discontinuities formed due to shock prop¬ 
agation are still not solved by FC methods and other 
additional numerical treatment must be applied for cor¬ 
rectly describe shock formation, e.g. intensive compu¬ 
tations by high order accurate weighted essentially non- 
oscillatory schemes (FC/WENO) 34 . Unlikely, the com¬ 
putational times increases by using those intensive com¬ 
putational techniques and the multiresolution analysis to 
detect discontinuities in the domain. 

The aim of the present work is to present a general¬ 
ization of the constitutive relations of nonlinear acous¬ 
tics including multiple relaxation processes in a non- 
convolutional formulation that allows the time-domain 
numerical solution by an explicit finite differences nu¬ 
merical scheme. Frequency power law attenuation based 
in relaxation have been applied in the same way than 
it has been applied to generalized Burgers equation 3 , 
Khokhlov-Zabolotskaya-Kuznetsov (KZK) model 3,4 and 
Westervelt equation 16 . The relaxation parameters have 
been fitted to both exact frequency power law attenua¬ 
tion and empirically measured attenuation of a variety 
of tissues that does not fit an exact power law. Two 
processes have been enough to model tissue attenuation 
with acceptable accuracy over a frequency range covering 
about 4 octaves, as it was demonstrated by Yang et al. 4 . 
A numerical technique based on artificial relaxation is in¬ 
cluded to control the non-negligible numerical dispersion 
of the FDTD method and improve stability when shock 
waves are present in the solution. The method includes 
backscattering and arbitrary propagation direction of fi¬ 
nite amplitude beams, and can be specially suitable in 
axisymmetric configurations where the computational re¬ 
sources for full 3D fc-space methods are prohibitive. 

The paper is organized as follows: in Sec. II the model 
equations that describes the problem are exposed, Sec. Ill 
describes the computational method presented in this 
work and in Sec. IV the method is validated comparing 
the numerical results with analytic solutions for linear, 
smooth and discontinuous nonlinear waves. In Sec. V so¬ 
lutions for frequency power law media are presented and 
compared with analytic and numerical solutions obtained 
by fc-space methods as benchmark case. Furthermore, an 
experimental test was presented where a focused beam 
propagating in castor oil has been modeled. Finally, a 
high intensity focused source in the high nonlinear regime 
focusing in soft-tissue is modeled. 


II. GENERALIZED NONLINEAR ACOUSTICS MODEL 
FOR MULTIPLE RELAXATION MEDIA 

A. Full-wave modeling 

The principles of mass and momentum conservation 
lead to the main constitutive relations for nonlinear 
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acoustic waves, which for a fluid can be expressed as 35 

Tt = “ v ' ( ' ,v) (1) 


and 

P + v • = -Vp + r;V 2 v + (c + |) V (V • v), 

( 2 ) 

where p is the total density field, v is the particle velocity 
vector, p is the pressure, and 77 and £ are the coefficients 
of shear and the bulk viscosity respectively. The acoustic 
waves described by this model exhibit viscous losses with 
quadratic power law dependence on frequency. In order 
to include a power law frequency dependence on the at¬ 
tenuation, a multiple relaxation model will be added into 
the time domain equations. 

The basic mechanism for energy loss in relaxing media 
is the appearance of a phase shift between the pressure 
and density fields. This behavior is commonly modeled 
as a time dependent relation at the fluid state equation, 
that for a fluid retaining the material nonlinear effects 
up to second order an be expressed as 35,36 : 


Vn = (c 2 — Cq)/cq, where c n is the sound speed in the 
high frequency limit associated to n-th order relaxation 
process, also known as the speed of sound in the “frozen” 
state 37 . In order to describe relaxation without the need 
of including a convolutional operator, we shall define a 
state variable S n for each process as 

Sn = —G n * p'. ( 6 ) 

T n 

Thus, using the convolutional property 

jft i G (t) * p'{t)) = * p'{t) = G(t) * the 

time derivative of the relaxation state variable obeys the 
following relation for the n-th order process: 


dS n / 
dt V 


1 i] 


n 0 —- 


' H(t) + 7 -^^-c <5(t)^ * p\ (7) 

7n / 


where 8{t) is the Dirac delta function. Using the Eq. ( 6 ) 
this relation becomes a simple ordinary differential equa¬ 
tion for each process as: 


dS n _ 1 p n cl , 

a, — ‘- > n + P 

dt T n T n 


( 8 ) 


p " c2 ° p ' + C %YX p2 + G[t ~ 


( 3 ) 


where p’ = p — po is the density perturbation over the 
stationary density po, B/A is the nonlinear parameter, 
Co is the small amplitude sound speed, and G(t) is the 
kernel associated with the relaxation mechanism. The 
first two terms of the right hand side of Eq. (3) de¬ 
scribe the instantaneous response of the medium, where 
the convolutional third term accounts for the “memory 
time” of the relaxing media. Thus, by choosing an ad¬ 
equate time function for the kernel G(t) the model can 
present an attenuation and dispersion response that fits 
the experimental data of the heterogeneous media. How¬ 
ever, the direct resolution of the constitutive relations 
Eqs. (1-3) in this integral form is a complex numerical 
task due to the convolutional operator. Thus, instead of 
describe G(t) with a specific time domain waveform, the 
response of the heterogeneous medium can be alterna¬ 
tively described by a sum of N relaxation processes with 
exponential time dependence as: 


Using again convolutional properties, we can substi¬ 
tute Eq. ( 8 ) into Eq. (4), and the relaxing nonlinear state 
Eq. (3) becomes: 


p = c 0 p + 


c\B_ 
Po 2^4 


N 


N 


p' 2 - Sn + 7inC o p' ■ ( 9 ) 


Moreover, if “frozen” sound speed for N mechanisms 

( N \ 

is defined as €^ = 00 ( 1 - 1 - Vn) > Eq. (9) leads to: 


P = clp' + C f^p' 2 -Y, S n- ( 10 ) 

Due to the smallness of the relaxation parameter, rj n , 
i.e. when weak dispersion is modeled, the sound speed in 
the high frequency limit reduces to 35 : 


N 


N 


Cqo — Cq 


i + E? ■ 


71=1 


(ID 


=!>*§(’ « 4 » 

with the n-th order relaxation kernel expressed as 

G n {t) = r) n cle^H(t), (5) 

where H(t) is the Heaviside piecewise function H(t < 
0) = 0, H(t > 0) = 1, r„ is the characteristic relaxation 
time and ?/„ the relaxation parameter for the n-th order 
process. This last dimensionless parameter controls the 
amount of attenuation and dispersion for each process as 


Note Eq. (10) for a mono-relaxing media is equivalent to 
that can be found in literature 36 


t 



— OO 





( 12 ) 


Thus, the constitutive relations to solve by means of 
the numerical method in the nonlinear regime are the 
continuity Eq. (1), the motion Eq. (2) and the second or¬ 
der fluid state relaxing Eq. (10), where the state variable 
S n obeys the relation Eq. ( 8 ) for the n-th order relax¬ 
ation process. Although the aim of this work is to model 
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biological media, the generalized formulation presented 
here can be used to describe the attenuation and hence 
the dispersion observed in other relaxing media, as the 
relaxation processes of oxygen and nitrogen molecules in 
air or the relaxation associated with boric acid and mag¬ 
nesium sulfate in seawater 37 . 


B. Small amplitude modeling 

If small amplitude perturbations are considered, an 
equivalent derivation of this model can be expressed for 
multiple relaxation media 25 . Thus, for an homogeneous 
inviscid relaxing fluid, the linearized continuity and mo¬ 
tion Eq. (1-2) reduces to 


Thus, the linearized governing Eq. (14, 18) for a relax¬ 
ing media are expressed in a pressure-velocity formula¬ 
tion and can be solved together with the coupled state 
evolution Eq. (19) by means of standard finite differences 
numerical techniques 25 . In this way, lossless linear acous¬ 
tics equations can be obtained by setting r] n = 0 or in the 
limit when the relaxation times r n —> oo. The relaxation 
behavior described by this linearized model is achieved 
too by the formulation described by 25 , where the relax¬ 
ation coefficients ry n and the relaxation variable S n are 
defined in a different, but analogous way. 

III. NUMERICAL SOLUTION BY FINITE-DIFFERENCE 
TIME-DOMAIN 


^ = -A)V-v (13) 

and 

<9v 

Po-gj- = -Vp; (14) 

and linearizing the fluid state Eq. (10) we obtain: 

p'=jr (p +£ ( is ) 

00 \ n=1 J 

These equations can be solved directly in this form, 
however, if expressed in pressure-velocity formulation the 
density field is no longer necessary and computational ef¬ 
fort can be reduced. Thereby, assuming a linear “instan¬ 
taneous” compressibility Koo = Poc^, and substituting 
Eq. (15) into Eq. (13) yields 


9p W dSn 
dt + ^ dt 

n —1 


-KooV • V. 


(16) 


Then, taking the time derivative of the state variable 
Eq. (8) we get 


dp 

m 


N 

£7 

1 ' 1 
71=1 


-s n +,/f 2 — 

, TYi 
71=1 


— K00V • V. 


(17) 


Finally, substituting again the linearized state Eq. (15) 
and arranging terms the linearized continuity equation 
leads to 


dp t ^ p n cl 


dt P ^ r n cL 


N 

£ 

71=1 


VnCl 

T n cL 


-] S rl = t 


(18) 

On the other hand, the state evolution equation can 
be expressed as a function of the acoustic pressure as 


dS n 

dt 



Pncl 

TnClo 



(19) 


In this section the numerical techniques for solving the 
complete set of equations (continuity Eq. (1), momentum 
Eq. (2), state Eq. (10) and the relaxation Eq. (8)) are 
presented. The numerical method is based on a second 
order FDTD method where multiple relaxation processes 
are included in order to: first, modeling physical atten¬ 
uation and dispersion at the frequencies of interests and 
second, correct the numerical dispersion and include ar¬ 
tificial attenuation to guarantee convergence in nonlinear 
regime. Moreover the inclusion of relaxation processes in 
the presented formulation require only one extra field per 
relaxation process and no memory buffer is needed. 


A. Discretization 


Cylindrical axisynnnetric x = (r, z) coordinate system 
is considered in this work, however, the method can be 
derived in other coordinate systems. As in the standard 
acoustic FDTD method 38 , the particle velocity fields are 
discretized staggered in time and space respect to the 
density and pressure fields. As shown in Fig. 1 uniform 
grid is considered, where r = zAr, z = jAz, t = mAt, 
with Ar and A z as the radial and axial spatial steps, and 
At is the temporal step. 

Centered finite differences operators are applied over 
the partial derivatives of the governing equations. Thus, 
spatial interpolation is needed over the off-center grid 
variables in order to fulfill the conservation principles 
over each discrete cell of the domain 39 . The r compo¬ 
nent of Eq. (2) is expressed in a cylindrical axisynnnetric 
system as 


dv r 

~dt 


1 dp 




p dr 


+ 


d 2 


dv r 

lfr~ Vz ' 

1 dv r 
dr 2 + r dr 

r , 1 \ (d 2 v r 
C+ 3V \~d^ 


dv r 
dz 

d 2 i v v r 

dz 2 r 2 
■ . 1 dv r + d 2 v z 
r dr drdz 


( 20 ) 



Each term of the above expression is approximated by 
centered finite differences evaluated at r = (i + |) • Ar, 
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FIG. 1. Spatial staggered discretization. The pressure (pY)j) 
and the n-th order relaxation process state fields (S^ij) are 
evaluated at same discrete location as the density (pYj). Par¬ 
ticle velocity fields are discretized staggered in both space and 
time respect to the density, pressure and the n-th order state 
fields. 


z = (j + \) • A z, t = (m + i) • At. This equation can 

be solved obtaining an update equation for 2 .. In 

the same way, the z component of the motion Eq. (2) is 
expressed as 


dv z 

~dt 


1 dp 
p dz 


+5 


+ - C+ nV 


dv z 

dr ~ Vz 

_1 dv z 

dr 2 r dr 

I„I ( a 2 ” 


d 2 


dv z 
dz 
d 2 v z 

I- 4 

dz 2 
- 1 dv 


3 ) \dzdr 


d 2 v z 
r dz ^ dz 2 


( 21 ) 


This equation is approximated by centered finite dif¬ 
ferences and evaluated at r = i ■ Ar, z = (j + d ) • Az, 
t = m ■ At. An update equation is obtained solving this 

equation for v z . . \ . Equation (1) in cylindrical axisyrn- 
*J+2 

metric coordinate system is expressed as 


dp 

dt 



dv z \ dp dp 

!h ) ~ Vr lfr~ Vz lh' 


( 22 ) 


Following the same procedure, each term of the above 
expression is approximated by centered finite differences 
and evaluated at r = i-Ar, z = j-Az, t = (m+|)-At, and 
the update equation is obtained solving this expression 
for p™ +1 . A leap-frog time marching is applied to solve 
Eq. (20-22) for each time step until the desired simulation 
time is reached. Finally, Eq. (8) is locally solved for to+ 1 
by and explicit fourth-order Runge-Kutta method and 
then Eq. (10) is used for update the pressure held. 


B. Boundary conditions 


nodes, allowing to prevent the singularity of the cylindri¬ 
cal coordinate system: due to the staggered grid, the 
only variable discretized at r = 0 is v r , and axisym- 
metric condition v r \ r —o = 0 is applied there. Further¬ 
more, to solve spatial differential operators at boundaries 
some “ghost” nodes must be created with the conditions: 
v r (—r) = —v r (r), v z (—r) = v z (r), p(-r) = p{r) and 
p(-r) = p(r). 

Perfectly matched layers (PML) 40 were placed in the 
limits of the domain (±z and +r) to avoid spurious reflec¬ 
tions from the limits of the integration domain. Inside the 
PML domains linearized acoustic equations were solved 
using the complex coordinate screeching formulation 41 . 
For a layer of 30 elements and a broadband incident 
wave with 1 MHz central frequency and non-normal in¬ 
cidence angle, these absorbent boundary conditions have 
reported a reflection coefficient of R = —55.2 dB. How¬ 
ever, the performance of the PML is amplitude depen¬ 
dent as long as the nonlinear terms are uncoupled to the 
PML domains. The amplitude dependence of the reflec¬ 
tion coefficient is shown in Fig. 2, where a PML of 25 
layers have reported reflection coefficients R < — 50 dB 
for waves in linear regime and highly nonlinear waves 
including shocks. 


C. Minimizing numerical dispersion 

The stability for the lossless linear FDTD algorithm 
follows the Courant-Friedrich-Levy (CFL) condition, 
that for uniform grid (A?’ = Az = Ah) the maximum 
duration of the time step is limited by At < Ah/cQ\fD 
where D is the number of dimensions (i.e. D = 2 in cylin¬ 
drical axisymmetric coordinate system). That condition 
essentially states that for a single time step information 
can not propagate in the numerical grid a distance longer 



The staggered grid is terminated on velocity nodes, so FIG. 2. Reflection coefficient of the perfectly matched layer 

the boundary conditions are applied on these external ( PML ) versus la y er thickness for different wave am P litude s- 
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FIG. 3. Normalized dispersion relation for of a FDTD lattice 
of Courant number S = 0.9 and anomalous dispersion relation 
for a mono relaxing acoustic media. The straight line Co = 
oj/k represents the reference nondispersive case. 


that one cell. However, if relaxation is included numer¬ 
ical instabilities have been observed when Tf /2tt < At. 
Due to this empirical relation, the maximum values for 
relaxation frequencies are limited too by the chosen spa¬ 
tial discretization by the simple relation f n < \/2N\f 0 , 
where /„ = 2tt /r n is the maximum relaxation frequency 
for all processes, N\ is the number of spatial samples 
per wavelength and /o the frequency of the propagating 
wave. 

On the other hand, nonlinear effects induce the pro¬ 
gressively growing of harmonics of the fundamental fre¬ 
quency of the initial wave. The diffusive viscous terms 
in Eq. (2), attenuates the small-amplitude high-spatial 
frequencies, damping the “node to node” numerical os¬ 
cillations and ensuring numerical stability in weakly non¬ 
linear regime. Thus, for a smooth solution the numerical 
algorithm shows consistency when Ah —> 0, so if stabil¬ 
ity is achieved by the CFL condition, the convergence 
is guaranteed. However, in strongly nonlinear regime, 
i.e. when sharp waveforms or even shocks are present 
in the solution, extra numerical techniques must be em¬ 
ployed to guarantee convergence. Artificial viscosity can 
be added when shock waves are present in the solution 
where a common implementation follows a fourth order 
spatial filtering 28,29 . Thus, the artificial attenuation re¬ 
trieved by this spatial operator is fourth power of fre¬ 
quency: the low frequency components of the solution 
remains quasi-undamped, while the higher spatial fre¬ 
quencies are strongly attenuated. In this way, the so¬ 
lution is smoothed and shock thickness depends on the 
artificial viscosity coefficient. 

However, the main drawback for finite difference meth¬ 
ods is numerical dispersion, where the analytic disper¬ 
sion relation can be expressed in ID as sin 2 = 


■jTsin 2 ^^^, with the Courant number S = CgAt/Ah. 
In this way, numerical dispersion reduces phase speed for 
high frequency components so traveling sharp solutions 
develop tail oscillations: the low wavenumbers travels 
fast and left behind high spatial frequencies. In nonlin¬ 
ear regime, is well-known that the combined effects of 
nonlinearity and strong dispersion can lead rich phenom¬ 
ena, e.g. beatings in the generated harmonics, pulsations 
on the vertex of a sawtooth wave or soliton formation in 
strong dispersive media 36 . In this way, the numerical dis¬ 
persion by discreteness of the FDTD methods couple to 
the physical nonlinearity can lead to a great variety of 
non-physical or even unstable solutions. 

In order to overcome those two limitations, i.e. the 
generation of harmonics over the discrete limit and the 
numerical dispersion, we propose the use of artificial 
relaxation. As Fig. 3 shows, physical relaxation pro¬ 
cesses introduces anomalous dispersion, i.e. the phase 
speed increases in the high frequency regime, opposite 
to the numerical (lattice) dispersion of the finite differ¬ 
ences scheme. Thus, by introducing a collection of re¬ 
laxation processes and choosing its adequate relaxation 
parameters the high frequency numerical dispersion can 
be compensated. As a consequence, introduction of those 
relaxation processes in the high frequency lead to the in¬ 
evitable inclusion of artificial attenuation. However, this 
numerical attenuation is then exploited to limit the grow¬ 
ing of higher harmonics in a similar way than artificial 
viscosity 29 . It is worth noting here that, due to the at¬ 
tenuation using artificial relaxation is, at maximum, only 
second power of frequency, the low frequency range of the 
solution is therefore also attenuated. Thus, the proposed 
method is restricted to lossy media. 

The adequate relaxation parameters that corrects the 
numerical dispersion have been found by multi-objective 
optimization techniques, where two cost function are pro¬ 
posed: one for dispersion and other for attenuation. In 
the first case, the error between the goal (ideal) disper¬ 
sion relation and the retrieved numerical dispersion cor¬ 
rected by multiple artificial relaxation and is evaluated 
in the high frequency regime. Finally, the second cost 
function is the error between the desired (ideal) attenu¬ 
ation and the artificial attenuation evaluated in the low 
frequency regime. 


IV. VALIDATION 
A. Single relaxation process 

A canonical case of a physical single relaxation process 
is presented. In order to correct numerical dispersion 
the parameters of two extra artificial relaxation processes 
have been found using the multi-objective genetic algo¬ 
rithms provided by the optimization toolbox in MAT- 
LAB R2014a v8.03. Linear propagation was considered 
and simulation parameters were set to typical values for 
water: cq = 1500 m/s, po = 1000 kg/m 3 , B/A = 5, 




7 


rj = 8.90 • 10~ 4 Pa-s. A single physical relaxation pro¬ 
cess was included, with a characteristic relaxation time 
of 7~i = 1/27t/o and fo = 2 MHz, and relaxation modu¬ 
lus of r]i = 0.0678 that leads to a frozen sound speed of 
Coo = 1550 m/s. In this case, the numerical parameters 
were set to Ar = Az = 1.87-10 -7 nr and At = 8.65-10” 11 
s. A plane wave traveling in +z direction was considered. 

Thus, the theoretical attenuation for the relaxation 
processes and including viscosity can be expressed as 


a(w) 


( t , f \ A iin u 2 t 2 

2 Po4 V 3 / 2cor n 1 + u 2 t 2 ’ 


(23) 


and the theoretical phase velocity can be predicted as 37 


Cp(uj) = c 0 



r] n uj 2 t l A 

2 1 + 


(24) 


In order to compute the attenuation and dispersion of 
the numerical method, simulated pressure was recorded 
at two locations Zq and z ±, and attenuation and phase ve¬ 
locity were estimated from the spectral components over 
the bandwidth of the input signal. The numerical atten¬ 
uation was calculated as 


ot{u>)f d 


ln(|P(q;,zi)/P(a;,zo)|) 
(zi - z 0 ) 


(25) 



FIG. 4. (Color online) Pareto front retrieved by the multi¬ 
objective genetic optimization. Square marker (A) is the so¬ 
lution those relaxation parameters minimizes the numerical 
dispersion. The best fit in the attenuation are the parameters 
that provided the solution marked by the triangle marker (B). 
A compromise between attenuation and dispersion is achieved 
at the individuals around the center of the Pareto front, as 
shows the sample marked by the circle (C). The inset shows 
the normalized dispersion relation retrieved by the individu¬ 
als (A) dotted line, (B) dashed-dotted line, (C) continuous 
line. Dashed black line shows the numerical dispersion re¬ 
lation of the FDTD method for a Courant number of 0.94 
without artificial relaxation. 


where P(w) is the Fourier transform of the measured 
pressure waveforms at points zq and z\. On the other 
hand, the phase velocity was computed as 


c p( w )/d — 


U‘{Z\~ Zq) 

arg (P(w, Zi)/P(uj, Zq)) ’ 


(26) 


where correct phase unwrapping is needed in the arg 
function. 

In this way, Fig. 4 shows the retrieved Pareto front of 
the optimization, where 3 different areas can be distin¬ 
guished. The first area, marked as (A) in Fig. 4, rep¬ 
resents individuals whose numerical dispersion is mini¬ 
mal but attenuation is not optimal. On the other hand, 
the individuals around area (B) represent a set of relax¬ 
ation parameters that provides the best agreement be¬ 
tween numerical and physical attenuation. A good com¬ 
promise between both situations can be obtained in the 
central area of the Pareto front (C), where retrieved at¬ 
tenuation and dispersion in the low frequency band shows 
good agreement with the physical, and the numerical dis¬ 
persion has been corrected over a wide frequency range. 
However, as can be seen in the inset of Fig. 4, the disper¬ 
sion relation retrieved by all the cases corrects the FDTD 
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FIG. 5. (Color online) Retrieved attenuation (top) and dis¬ 
persion (bottom) by the inclusion of artificial relaxation for 
the individuals A (dotted), B (dotted-dashed) and C (dashed) 
of the individuals marked in Fig. 4. Continuous lines repre¬ 
sent the frequency range included by the optimization, while 
dashed lines shows the not optimized frequency range. As 
a consequence of correcting dispersion, attenuation increases 
in the high frequency range, as shown in the top subplot. 
This extra attenuation is used as artificial attenuation for 
numerical-nonlinear stability. 





















lossless numerical dispersion relation. The phase speed 
of those tree individuals is shown in Fig. 5 b), where it 
can be seen that in the frequency range selected for the 
optimization the numerical phase speed is corrected for 
all the individuals, where the best fit is obtained for in¬ 
dividuals in the Pareto front area A. On the other hand, 
the inclusion of artificial relaxation leads to an increas¬ 
ing of the attenuation in the high frequency range, as is 
shown in Fig. 5 a). In this way, as the phase speed error 
is reduced the effect of artificial attenuation increases. 
Although this increasing can be seen as a non-desired 
counterpart, the appearance of this attenuation is use¬ 
ful in order to control the harmonic growing in nonlin¬ 
ear regime in the same way as artificial viscosity spatial 
operators 29 . 


B. Nonlinear steady solution for single relaxation process 


In order to validate the method in the nonlinear regime 
a full-wave simulation was developed in a mono-relaxing 
media using above parameters. Thus, the analytical (in¬ 
verted) solution for the steady solution with p = —p o for 
r = — oo, p = po for t = oo and p = 0 for r = 0, for the 
retarded time r = t — z/cq reads 22 


D -1 


= T n In 


(1+p/po) 
(l^p/p 0 ) D+1 


(27) 


where D = rj n poCQ/2/3po measures the ratio of relaxation 
effects to nonlinear effects. For D > 1, where no shock is 
present, the solution retrieved by FDTD algorithm shows 
good agreement with analytical and no artificial attenu¬ 
ation is needed. However, for D < la discontinuity is 
present in the solution and convergence is only possible 
with the inclusion of extra numerical techniques. 

Thus, Fig. 6 (a-c) shows the analytical and numerical 
solutions including artificial relaxation and artificial vis¬ 
cosity, where excellent agreement is achieved in all cases. 
In the case of artificial viscosity, higher harmonics are 
strongly attenuated and by reducing grid step conver¬ 
gence can be achieved. Due to artificial viscosity operator 
is essentially a low pass spatial filter, a smoothed version 
of the shock is achieved. However, the phase speed of the 
higher spatial frequencies present in the shock is modified 
due to numerical dispersion, and for D « 1 (Fig. 6 (b, 
c)) oscillations appears in the tail of the discontinuity, 
leading to the appearance of non-physical solutions. 

On the other hand, the proposed method of artificial 
attenuation by relaxation also limits the harmonic grow¬ 
ing so a smooth version of the shock appears. More¬ 
over, artificial relaxation also corrects phase velocity so 
all the spatial frequencies travels at same speed and no 
oscillatory tail appears. The case of D = 0.01 is shown 
in Fig. 6 (c), where nonlinear effects strongly dominates 
over attenuation. In this case, tail oscillations provided 
by artificial viscosity increases in amplitude. In contrast, 
by including artificial relaxation a smoothed version of 
the shock is captured and accuracy is maintained. 



r/Dt r 




FIG. 6. (Color online) (a) Analytical (thick gray line), nu¬ 
merical using artificial relaxation with corrected dispersion 
(continuous line) and using artificial viscosity 29 (dotted line) 
for the nonlinear steady state solution for D = 1. (b) Non¬ 
linear steady state solution for D = 0.1. Inset shows detailed 
shock numerical solution for artificial relaxation (x markers) 
and artificial viscosity (+ markers), (c) Nonlinear steady state 
solution for D=0.01. 


V. RESULTS 

A. Frequency power law attenuation 

Using methodology described above, the optimal re¬ 
laxation parameters were obtained in order to fit the 
multiple-relaxation numerical attenuation to frequency 
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FIG. 7. Attenuation retrieved by the numerical algorithm 
(markers) and target frequency power law attenuation (gray 
lines). By using the optimization algorithm the relaxation 
times and modulus were optimized for minimize the relative 
error between the target power laws of 7 = [1, 1.3, 1.6, 2] and 
the attenuation retrieved. 


FIG. 8 . Phase speed retrieved by the numerical algorithm 
(markers) and target frequency power law attenuation (gray 
lines) for 7 = [1,1.3,1.6, 2], 


B. Fitting attenuation for tissue experimental data 


power law attenuation in the form 

a(oj) = crow 7 , (28) 


where 7 is the power law exponent and ao the power law 
coefficient in Np (rad/s) 7 m _1 . Moreover, the numerical 
dispersion was corrected by means of artificial relaxation 
in order to fit the corresponding frequency power law 
dispersion, where its analytical form satisfying causality 
can be expressed as 42 


Cp 


1 

(w) 


1 

c(w 0 ) 


+ a 0 tan (|w | 7 1 



This expression is valid for 0 < 7 < 3 with 7 / 1, 
and an alternate equation can be found 42 in the limit for 
7 = 1. Here, simulation parameters were Co = 1500 m/s, 
po = 1000 kg/m 3 , B/A = 5, 77 = 8.90 • HT 4 Pa-s, / 0 = 1 
MHz, Ar = Az = 1.3 • 10 -5 m, At = 5.4 • 10 ~ 9 s; that 
leads to 26 elements per wavelength and a CFL number 
of 0.9. Only two independent relaxation processes were 
employed in this section to obtain the target frequency 
power laws. 

Following the above procedure, the relaxation times r n 
and relaxation modulus r) n were optimized for different 
frequency power laws covering the range of that observed 
in tissues 7 = [1,1.3,1.6, 2]. The attenuation coefficient 
ao was chosen for each power law to present an atten¬ 
uation a = 1 dB/cm/MHz 7 . The fitting was developed 
over the typical frequency range for medical ultrasound 
applications, i. e. 1 to 20 MHz for both attenuation and 
phase speed. The results for the attenuation and phase 
speed curves are plotted in Fig. 7 and Fig. 8 , where the 
theoretical and the numerical predictions agree over the 
frequency range used for the fitting. 


Although a frequency power law dependence can de¬ 
scribe the ultrasound attenuation over a finite frequency 
range, the attenuation data of some particular examples 
shows variation of the exponent over the entire frequency 
range 43 . Thus, Fig. 9 shows experimental attenuation 



FIG. 9. Experimental attenuation data for some tissues 
adapted from Ref. Hill, Bamber, and ter Haar 43 (lines), and 
obtained by the numerical method (markers) by fitting the 
parameters of 2 relaxation processes. The numbers above the 
curves show the exponent of the frequency power law 7 for 
each frequency region (i.e. the slope of the curve). 
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TABLE I. Error of the optimized attenuation response rela¬ 
tive to the experimental data for N total relaxation processes. 


Tissue 

Power law 

N = 1 

AT = 2 N = 3 

N = 4 


(local slope) 

£(%) 

£(%) 

£(%) 

£(%) 

Skin 

f os 

6.67 

0.167 

0.136 

0.120 

Liver 

f 1 - 2 

7.62 

0.517 

0.404 

0.165 

Blood 

z 1 - 4 , f 1 

8.34 

0.349 

0.330 

0.310 

Breast 

/°' 9 , f 1 - 2 

5.20 

0.216 

0.209 

0.205 

Skull bone /°' 9 , / 21 , /°' 5 

10.60 

10.54 

8.628 

5.189 


data curves for some tissues where the local slope of the 
power law changes over the measured frequency range. 
This behavior can be modeled by a sum of relaxation 
processes by optimizing the relaxation parameters as de¬ 
scribed above. Thus, the results show that most tissues 
with locally variable 7 can be fitted by only a pair of 
relaxation processes, as the same way that for constant- 
slope frequency power law attenuation 3 . 

In this way, Table I shows the error of the nu¬ 
merical attenuation relative to the experimental data. 
The percent relative error was computed as e = 
T^Ti ffi df, where a e (f) is the experimental 

attenuation data, f± and /2 define the frequency range of 
the measurement. 

As expected, the goodness of fit grows as the number 
of relaxation processes included increases. However, only 
two processes are enough to obtain relative errors below 
1% for tissues with 7 < 2. In the case of tissues where a 
local value of 7 > 2 has been observed, the fitting pro¬ 
cedure fails, like in the skull bone in the 2 MHz range 43 . 
The maximum slope achieved by single relaxation and 
thermo-viscous losses is 7 = 2 for any frequency, so a tis¬ 
sue showing that slope cannot be accurately modeled in 
this frequency region with the method proposed in this 
work. From another point of view, Eq. (29) states that 
frequency power law medium with 2 < 7 < 3 presents 
standard dispersion 42 , opposite to anomalous dispersion 
for media falling in the range 0 < 7 < 2. Therefore, the 
dispersion relation of media with 2 < 7 < 3 cannot be 
modeled by a sum of relaxation processes as long relax¬ 
ation includes only anomalous dispersion. 


TABLE II. Variation of sound speed (Ac) observed numeri¬ 
cally for the modeled tissues by means of two relaxation pro¬ 
cesses and analytical using the Kramers-Kronig relations. 


Tissue 

Numerical Ac Analytical Ac 


(m/s) 

(m/s) 

Skull bone 

80.737 

70.720 

Skin 

10.148 

2.460 

Breast 

2.323 

2.455 

Liver 

3.118 

2.339 

Blood 

0.865 

0.907 


Using Kramers-Kronig relations 44 , the variations of 
sound speed Ac can be predicted by the frequency depen¬ 
dent attenuation. Table II shows the variation of sound 
speed observed in the numerical solution over the fitted 
frequency range. The magnitude of these variations are of 
the order of magnitude of those measured experimentally 
in this frequency range, and the frequency dependence 
observed for the variation is roughly linear as observed 
in real tissue 43 . As expected from the relations between 
dispersion and absorption 44 , the magnitude of the varia¬ 
tion in sound speed increases as the total variation of the 
absorption increases for a given frequency range. 


C. Nonlinear one-dimensional propagation in tissue-like 
media 

1. Non-dispersive media 

In order to study the convergence of the numerical cal¬ 
culations to an analytical solution of the model in the 
nonlinear regime, a medium with frequency squared de¬ 
pendence attenuation is implemented using the adequate 
relaxation times and relaxation modulus as explained 
above. The numerical solution is compared with the 
analytical solution for a plane wave traveling through a 
thermo-viscous fluid proposed by Mendousse 37 : 

P_ _ r S”=i (~l)" +1 dn ( 7 ) ne~ n2a / T sin {nut') 

Po Io ( 7 ) + 2 J2n=i (- 1 )"4 ( 7 ) ne _n2<T / r cos {nut') ’ 

(30) 

where T is the Gol’dberg number, defined for power 
law media as T = x a /x s , with absorption length x a = 
1 /aow 7 , shock formation distance x s = l/fiek and nor¬ 
malized distance a = x/x s \ with the parameter of non¬ 
linearity (5 = 1 + B/2A and the acoustic Mach number 
£ = v/co, with v is the source particle velocity and k the 
wavenumber. 

Figure 10 (top) presents the simulated waveforms at 
a = 1 and a = 3. The wave steepening due to nonlinear 
processes in the absence of dispersion are well resolved by 
the numerical method presented here. In order to study 
the accuracy of the algorithm, the amplitude of the first 
ten harmonics has been extracted for numerical and an¬ 
alytic solutions and plotted versus a in Fig. 10 (bottom). 
The observed relative error of the computational method 
decreases due to grid coarsening by a square law (i.e. 
the numerical scheme is second order accuracy). The 
magnitude of the error mainly depends on the number 
of elements per wavelength but, due to not ideal disper¬ 
sion, also on the traveled propagated distance. Including 
the correction of dispersion by artificial relaxation, for 
a path length of 100 A, a grid of 26 elements per wave¬ 
length was needed to obtain a relative error below 1 % 
for the third harmonic. Obviously, the relative error of 
the first and second harmonics will be always lower, i.e. 
the fundamental component error was 0.072 %. 
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FIG. 10. (Top) Waveforms at a = 1 and a = 3 for thermo- 
viscous attenuation (7 = 2 ). Mendousse analytical solution 
(black line), fc-space (gray line) and FDTD numerical solution 
(markers). (Bottom) Spatial distribution of the first ten har¬ 
monics for Mendousse analytical solution (black line), fc-space 
(gray line) and FDTD numerical solution (markers). 


FIG. 11. (Top) Waveforms at a — 1 and a = 3 in a tissue-like 
media with frequency power law (7 = 1 . 6 ). fc-space (gray line) 
and FDTD numerical solution (markers). (Bottom) Spatial 
distribution of the first ten harmonics for fc-space (gray line) 
and FDTD numerical solution (markers). 


In addition, the solution was compared also to the 
obtained by a fc-space method applied to the constitu¬ 
tive relations, i.e. the fc-wave algorithm 32 . This method 
was selected due to the low numerical dispersion and the 
possibility of including frequency power law attenuation. 
The result of both computational methods and Eq. (30) 
agree over all the spectral components analyzed, showing 
convergence to the analytic solution. 


2. Dispersive media 

In the case of frequency power law attenuation me¬ 
dia with 7 = ( 1 , 2 ) no general analytic solution exist in 
nonlinear regime for monochromatic progressive waves. 
Thus, in order to study convergence in this regime, the 
proposed FDTD solution was compared with the solution 


obtained by fc-space methods 32 . By using same physi¬ 
cal and grid parameters in both methods, the solutions 
agrees for different power laws. Thus, Fig. 11 (top) shows 
the good agreement for the waveforms measured at a = 1 
and cr = 3 obtained for 7 = 1.6, where the characteris¬ 
tic asymmetry effect of media with anomalous dispersion 
(e.g. relaxing, boundary layer effects ) 22 is observed: the 
shock front after the rarefaction phase is followed by a 
rounded positive compression profile. The spatial distri¬ 
bution for each harmonic is shown in Fig. 11 (bottom), 
where it is observed that the proposed FDTD solution 
with optimized attenuation and dispersion converges to 
the obtained by pseudo-spectral methods up to ten har¬ 
monics. As in the case of frequency squared media, grid 
refinement numerical tests have reported a second order 
accuracy of the FDTD method in nonlinear regime. 
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D. Nonlinear propagation in tissue-like media including 
diffraction 

1. Experimental validation 

An experiment was designed to test the validity of the 
algorithm for intense beams in frequency power law at¬ 
tenuation media. The source was formed by a plane 
single element piezoceramic crystal (PZ 26, Ferroperm 
Piezoceramics, Denmark) mounted in a custom designed 
steel housing and a polymethyl methacrylate (PMMA) 
focusing lens with aperture A = 50 mm and radius of 
curvature R = 50 mm. The source was driven with 
a sinusoidal pulse burst of frequency /o = 1.112 MHz 
and n = 50 cycles using a function generator (14 bits, 
100 MS/s, model PXI5412, National Instruments) and 
a linear RF amplifier (ENI 1040L, 400W, 55dB, ENI, 
Rochester, NY). The pressure waveforms were acquired 
with a HNR 500 pm needle PVDF hydrophone (Onda 
Corp, CA), and a digitizer (64 MS/s, model PXI5620, 
National Instruments) was used. A three-axis micropo¬ 
sitioning system (OWIS GmbH, Germany) was used to 
move the hydrophone in three orthogonal directions with 
an accuracy of 10 /im. The amplitude frequency response 
of the hydrophone was compensated in post-processing 
but not in phase due to the absence of phase calibration 
for this equipment. 

The source was completely immersed in a castor oil 
tank (350 x 350 x 350 mm). We select this frequency 
power law attenuation media due to the low variability 
of its acoustic properties along existent literature 30,45 . 
Using a sound speed inside the bulk of the lens Ci = 
2711 m/s, and a sound speed of the castor oil of Co = 
1480 m/s (at 26 C room temperature), the effective lens 
geometrical focal is estimated as F = R/( 1 — co/cj) = 
110.1 mm, leading to a linear lossless gain of G = 13.4. 

On the other hand, a nonlinear simulation including 
diffraction and frequency power law attenuation with 
same parameters was carried out in a workstation (20 
cores Intel Xeon E5-2680 CPU, 2.8GHz with 256 GB 
RAM). The boundary conditions were implemented for 
a spherical focused ultrasound source. The castor oil pa¬ 
rameters at 26 C room temperature 45 , were Co = 1480 
m/s, po = 961 kg/m 3 , a = 0.4 dB/cm/MHz 7 , 7 = 1.69, 
B/A = 12.0. The grid parameters were Ar = Az = 29.6 
/jrn and At = 13.6 ns, leading to a CFL number S = 0.95 
and N\ = 50 elements per wavelength at fundamental 
frequency, note that this grid leads to N\ = 16 for third 
harmonic. 

The balance between nonlinear effects and power law 
attenuation can be estimated by using the Gol’dberg ra¬ 
tio, r = x a /x s where x s is the shock formation dis¬ 
tance and x a the media attenuation characteristic length. 
Thus, the amplitude of the source were selected to obtain 
a Gol’dberg ratio of T = 0.25 in order to let the frequency 
power law attenuation effects slightly dominate over non¬ 
linear effects. On the other hand, the ratio between 
diffraction effects and nonlinear effects can be described 



z/F 

FIG. 12. Spatial distribution of the fundamental (o), second 
(+) and third (x) obtained by the numerical (continuous line) 
and experimental methods (markers) for a focused transducer 
immersed in castor oil. 

by the so called Khokhlov number 35 as Nk = x s /xd, 
where Xd = ka 2 / 2 is the diffraction length and a the 
source radius. For the proposed test a Khokhlov number 
of N = 0.5 was selected to let the nonlinear effect slightly 
dominate over diffraction effects. The selected excitation 
pressure amplitude was po = 87.7 kPa. 

The results are summarized in Fig. 12, where axial 
pressure distribution for the fundamental, second and 
third harmonic are presented. A good agreement is found 
between simulations and the experimental test. Only far 
to the focal point the amplitude there exist differences be¬ 
tween computations and experiments, that can be caused 
by nonuniform vibration of the source 46 , boundary effects 
of the PMMA lens, or miss-alignment of the source axis 
and micro-positioning system orthogonal directions along 
the 100 mm axial measurement. 

The maximum amplitudes of the first harmonic were 
p e i = 0.6539 MPa for the experiment and p n \ = 0.6576 
for the numeric. The second harmonic peak pressure was 
p e 2 = 78.368 kPa and p n 2 = 76.272 kPa, and the third 
harmonic peak pressure p e 3 = 10.146 and p n 3 = 10.252 
kPa for the experimental and numeric respectively. The 
relative errors between numerical and experimental re¬ 
sults are 0.56 % for the fundamental frequency, 2.67 % 
for the second harmonic and 1.04 % for the third. No 
error estimation was done for the peak pressure due to 
the absence of a phase calibration of the hydrophone. 


2. Highly focused beam 

In order to test the algorithm in the very high nonlinear 
regime with realistic tissue parameters a focused bowl of 
geometrical focal F = 50 mm and aperture A = 50 mm, 
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FIG. 13. Axial spatial distribution of the peak compression 
(p + ) and minimum rarefaction ( p ~) pressure for a focused 
beam propagating through a liver tissue layer. (Top) weakly 
nonlinear propagation (F = 0.16) and (bottom) strong non¬ 
linear effects (r = 0.79). The insets show the waveforms 
recorded at the geometrical focal. 


driven at /o = 1 MHz was numerically tested. These 
parameters leads to a source gain G = 26.5 and a /- 
number= 1 , showing that source is beyond the paraxial 
limit. The media consist in two layers. The first layer, 
where the source was located, was water at 20 C with 
parameters Ci = 1482 m/s, p\ = 1000 kg/m 3 , B/A\ = 5, 
a\ = 2.17 x 10 -3 dB/cm MHz 7 , 71 = 2. At a distance 
Zi = 0.7 /F = 35 mm, a layer of human liver tissue was 
placed, therefore the focal spot is located inside tissue 
at a depth of 15 mm. Liver tissue parameters 43 were 
c 2 = 1597, p 2 = 1050 kg/m 3 , B/A 2 = 7.9, a 2 = -75 
dB/cm MHz 7 , y 2 = 1.5. In this case, the numerical grid 
parameters were Ar = Az = 29.6 /im and At = 11.9 ns, 
leading to a CFL number S = 0.9 and N\ = 50 elements 
per wavelength at fundamental frequency ( N\ = 16 for 
third harmonic). 

Figure 13 shows the spatial peak pressure profiles for 
different excitation amplitudes. In Fig. 13 (top) the pres¬ 
sure of the source was po = 0.18 MPa, while in the 
bottom figure was increased to po = 0.94 MPa. Thus, 
for the selected parameters, Fig. 13 (top) present re¬ 
sults for r = 0.16 and Nk = 0.62, so the attenua¬ 
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FIG. 14. (Color online) Spatial distribution of the peak com¬ 
pression ( p + ) and minimum rarefaction ( p ~) pressure for a fo¬ 
cused beam propagating through a liver tissue layer, (bound¬ 
ary marked with dashed line) for T = 0.79. Colorbars are in 
|p|/po units. 


tion effects dominates over nonlinear effects and non¬ 
linearity slightly dominates over diffraction effects. In 
this way, low asymmetry is observed between the pos¬ 
itive compression peak, p + , and the minimum rarefac¬ 
tion pressure distribution, p~. The calculated wave¬ 
form at z = F, shown in the inset of Fig. 13 (top), is 
weakly distorted. However, there exist differences be¬ 
tween its normalized peak amplitude p + /po = 25.7 and 
P~/Po = 21.05, and the source characteristic linear gain, 
G = 26.5. They are caused, in one hand by the atten¬ 
uation effects, where the value of the lossy linear gain 
observed was G a = p + /po = P~/po = 23.1 measured at 
z = F, i.e. the amplitude at the focal was reduced to 
87.2% of the lossless amplitude. On the other hand, the 
differences due to the asymmetry between compression 
and rarefaction cycles are caused by the combined effect 
of nonlinearity and focusing. 

If source amplitude is increased to po = 0.94 MPa, 
as shown in Fig. 13 (bottom), the ratio between attenua¬ 
tion and nonlinear effects is increased to T = 0.79. In this 
regime, nonlinear effects are almost of the same order of 
attenuation effects. On the other hand, increasing source 
amplitude while keeping same transducer parameters im¬ 
plies also the Khokhlov number changes to Nk = 0.12, so 
the nonlinearity clearly dominates over diffraction effects. 
In this regime, highly asymmetric pressure distribution is 
observed, where the values at focal point are p + = 49.05 
MPa and p~ = —14.91 MPa. 

Other typical nonlinear phenomena characteristic of 
high intensity focused sources can also be observed: for¬ 
mation of sharp shock front and its corresponding har- 
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monic generation, or, as Fig. 14 shows, narrowing of the 
beam for p + and broadening for p~ pressure distribu¬ 
tions. In addition, nonlinear focal shift, i.e. displacement 
of the peak pressure relative to the position of the linear 
peak pressure can be also predicted for tissue propaga¬ 
tion. In the case of T = 0.79 it was observed a nonlinear 
focal shift A F + = +1.05 mm and A+ _ = —1.03 mm for 
the p + and p~ pressure distribution respectively. 


VI. SUMMARY 

A general model based on the full constitutive rela¬ 
tions of nonlinear acoustics in relaxing media have been 
presented in a time-domain formulation which does not 
require convolutional operators. A numerical solution 
by means of finite-differences in time-domain have been 
obtained, showing that the theoretical attenuation and 
dispersion due to relaxation processes can be achieved 
by the numerical method with accuracy. These results 
can be also used to model typical relaxation process of 
other relaxing media (e. g. the processes observed in air, 
associated to the molecules of oxygen and nitrogen, or in 
seawater, associated to the relaxation of boric acid and 
magnesium sulfate). 

Moreover, a method for modeling frequency power law 
attenuation by means of multiple relaxation has been im¬ 
plemented in the constitutive relations. The proposed 
method can describe local variations of the exponent of 
the frequency power law, so an arbitrary attenuation 
curve in the range 0 < 7 < 2 can be modeled by means 
of the proper optimization of the relaxation coefficients. 
This feature of the presented method is an advantage 
when compared with most fractional derivatives meth¬ 
ods, where the attenuation follows an exact but unique 
frequency power law over the entire frequency range. A 
broad range of human tissues have been modeled and 
the goodness of the fit using from two to four relaxation 
processes has been discussed. 

Furthermore, a computational technique that exploits 
the anomalous dispersion of the relaxation processes is 
employed to mitigate the numerical dispersion of the 
finite-differences scheme. Thus, while phase speed is 
corrected by including artificial relaxation processes, its 
corresponding artificial attenuation is used to improve 
stability in the nonlinear regime. In this way, smooth 
and stable versions of shock waves have been obtained 
and compared with its analytic solution. Furthermore, 
the validity of the algorithm including diffraction have 
been tested with experimental measurements of a focused 
beam in castor oil. 

Due to the model is developed from the constitutive 
relations for nonlinear acoustics, most wave phenomena 
is captured. As a difference from the one-way models 
the proposed model implicitly includes multiple wave di¬ 
rection, and, due to the Lagrangian density of acoustic 
energy is implicitly included in the computation, multiple 
scattering and strong resonance effects can be accurately 


described. Moreover, unlike KZK and other parabolic 
approximations, the proposed model captures the diffrac¬ 
tion exactly, so for simulation of acoustic beams the field 
is not approximated only to the beam axis, but also in 
the near field and far to the beam axis and thus high 
focused devices can be simulated. 

The code has shown to be particularly appropriate 
if the problem to simulate presents axisymmetry, be¬ 
cause the constitutive relations for nonlinear acoustics 
are solved in a computational 2D domain, while standard 
k -space methods need to employ full 3D domains due to 
the poor convergence of the Fourier series at discontinu¬ 
ities (r = 0). This is the case, for example, of the focused 
ultrasound transducer simulated in Sec. VD2, where a 
full 3D solution will require huge computational resources 
and calculation times. Finally, due to the particle ve¬ 
locity vector and the acoustic density fields are solved 
implicitly by the code, this information can be used to 
estimate other relevant magnitudes as the full nonlinear 
intensity vector, the nonlinear acoustic radiation forces 
in these absorbing media or the acoustic streaming gen¬ 
erated in frequency power law attenuation fluids. 
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